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ABSTRACT 

The rotation states of kilometer sized near earth asteroids are known to be affected by the YORP 
effect. In a related effect, Binary YORP (BYORP) the orbital properties of a binary asteroid evolves 
under a radiation effect mostly acting on a tidally locked secondary. The BYORP effect can alter 
the orbital elements in ~ I0 4-5 years for a D p — 2 km primary with a D s — 0.4 km secondary at 
1 AU. It can either separate the binary components or cause them to collide. In this paper we devise 
a simple approach to calculate the YORP effect on asteroids and BYORP effect on binaries including 
Ji effects due to primary oblateness and the sun. We apply this to asteroids with known shapes as 
well as a set of randomly generated bodies with various degrees of smoothness. We find a strong 
correlation between the strengths of an asteroid's YORP and BYORP effects. Therefore, a statistical 
knowledge on one, could be used to estimate the effect of the other. We show that the action of 
BYORP preferentially shrinks rather than expands the binary orbit and that YORP preferentially 
slows down asteroids. This conclusion holds for the two extremes of thermal conductivities studied 
in this work and assuming the asteroid reaches a stable point, but may break down for moderate 
thermal conductivity. The YORP and BYORP effects are shown to be smaller than what could be 
naively expected due to near cancellation of the effects on small scales. Taking this near cancellation 
into account, a simple order of magnitude estimate of the YORP and BYORP effects as function of 
the sizes and smoothness of the bodies is calculated. Finally, we provide a simple proof showing that 
there is no secular effect due to absorption of radiation in BYORP. 



1. INTRODUCTION 

When NEA orbit around the sun they are constantly 
subjected to the sun's radiation. In equilibrium, the to- 
tal energy absorbed by the asteroid must be re-emitted. 
Yet, asymmetry in the asteroid's geometry results in a 
resid ual force that te nds to change the asteroid's mo- 
tion (|Rubincaml I2000D . This residual force, called the 
YORP effect, can significantly change the spin rate of 
a kilometer sized asteroids at a distance of I AU from 
the sun in I0 5 - 10 6 years. The YORP ef fect has been 
successfully measured for several NEA dTavlor et al.l 
[20071: iKaasalainen et all 120071; iDurech et all 120081) . The 
high abundance of fast rotators in the NEA family 
(|Pravec et all 120081 ) might be explained by the YORP 
effect. A rubble pile NEA might undergo fission if the 
YORP effect accelerates it beyond it's rotational breakup 
velocity (| Vokrouhlickv fc Capekl [20021 ); this mechanism 
might explain the formation of binary NEA. 

The differential acceleration between the two compo- 
nents of the binary is mostly due to the acceleration of 
the secondary due to its larger surface to volume ratio. 
A coherent effect also requires at least one of the com- 
ponents to be tidally locked which is more common for 
the secondary. The net radiation force acting on the 
secondary will produce a torque relative to the primary. 
This Binary YORP (BYORP) effect, first suggested by 

ICuk fc Burnsl (|2005D . evolves the orbit of the binary on 
fairly fast timescales (~ 10 5 years for a D p = 2 km 
primary with a D s — 0.4 km at 1 AU separated by 
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a = l.bDp). The spin rate and the obliquity of the as- 
teroid evolve due to the YORP effect. Similarly, the 
semi-major axis of the binary and its inclination relative 
to the orbital plane around the sun evolve by the BY- 
ORP effect. In addition the BYORP effect changes the 
eccentricity vector of the binary. We show that there 
are preferred end states for both effects. YORP tends to 
slow the spin rates, while BYORP tends to shrink the 
semi-major axis in the binary case. 

In §2 we introduce our model method and assumptions. 
In §3 we show that neighboring areas on the surface of the 
asteroid tend to have opposing effects. We provide order 
of magnitude estimates of the YORP effect as function 
of the sizes of the body and its smoothness. A method 
of precise calculation, as well as counting for thermal lag 
due to finite thermal inertia is shown in §4 and the results 
are discussed in §5. Analogous results for the BYORP 
effect are derived in §6. 

2. STRUCTURE MODELS AND COORDINATE SYSTEMS 

2.1. Modeling Method 

We model the asteroids by means of tessellation, where 
the asteroid's surface is described by a set of triangles. 
We neglect the effect of shadows cast by one facet over 
another and we assume that the emission from each point 
on the surface of the asteroid is isotropic. Since the or- 
bital period of the asteroid around the sun is much longer 
than its rotation time, it is assumed that there are no res- 
onances between the asteroid's orbit around the sun and 
the asteroid's revolution around itself. Constant density 
is assumed, although possible effects of density non uni- 
formity will be briefly addressed. 

2.2. Coordinate Systems 
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Two coordinate systems will be used throughout this 
paper: 

1. An inertial frame with axes labeled x,y,z. The 
z axis is perpendicular to the orbital plane of the 
asteroid around the sun. The x axis is chosen to 
coincide with the projection of the asteroid's spin 
vector on the xy plane. The origin is chosen to 
be the sun. This system will be referred to as the 
inertial system. 

2. The principle axes of the asteroid, labeled x',y',z', 
where z' is parallel to the spin vector. This system 
will be referred to as the asteroid system. 

3. SCALING OF THE RADIATIVE TORQUE 

In this section we discuss the net acceleration and to- 
tal torque that arises from the emission of radiation by a 
body of arbitrary shape. From symmetry, spherical bod- 
ies do not exhibit torque or acceleration, and we derive 
here the general scaling of the torque and acceleration as 
function of the roughness of the body, or its deviations 
from sphericity. We treat only the effect that arises from 
emission of radiation since iRubincam fc Paddackl (|20i0f ) 
have shown that there is zero secular change due to ab- 
sorption of radiation. 

We assume that the orbital time around the sun and 
its rotational period are non commensurate, and we can 
therefore perform the time average by averaging over the 
orbit and the spin sequentially. We find it more con- 
venient to first fix the angle of rotation of the asteroid 
around its axis and average the torque applied to the as- 
teroid during an orbit around the sun, and then average 
the torque as the asteroid revolves around its axis. 

In order to calculate the effect of YORP, the complete 
structure of the asteroid needs to be known. In this sec- 
tion we consider the scaling relations of YORP and pro- 
vide order of magnitude estimate of the effect. 

The change in the spin rate, s, of a homogeneous as- 
teroid scales according to: 



~p~R? 



(1) 



where $ is the solar radiation momentum flux given by 



$ = 



L, 



47rcd 2 Vl - e 2 



(2) 



Here L & is the solar luminosity, e is the eccentricity, c is 
the speed of light, d is the orbital semi-major axis, p is 
the density, and R is the length scale of the asteroid. The 
eccentricity dependence arises from averaging the torque 
over a heliocentric orbit. 

We construct simple models to account for the asym- 
metry of the asteroid. We choose n points randomly dis- 
tributed on the unit sphere and connect them to create 
a tessellation of small triangles that encloses the asteroi d 
(based on the Quickhull algorithm IBarber et al.l {1996), 
which produces In — 4 triangular facets for a given n) . 
This method of construction eliminates shadowing of one 
facet over another. For this body, we now calculate the 
radiation effects, and their scaling with n or with the de- 
viation of the shape of the body from a sphere. In the 
estimates below, we assume n~^> 1. 



We define the deviation of the asteroid from a sphere 
with the same volume by comparing the normalized dif- 
ference in their surface areas. 



d, = 



S - 47T7- 2 
47T?" 2 



(3) 



where 47rr 3 /3 is the volume of the asteroid and S is the 
surface area of the asteroid. This definition of spherical 
deviation will be shown to be d s oc n~ 1 . In order to 
simplify the derivation we will assume that all of the 
facets are equilateral and that the center of mass (CM, 
hereafter) is located at the origin. If we connect each one 
of the facets to the CM and thereby create a tetrahedron, 
the relation between 9, which is the vertex angle of a 
tetrahedron face, and £1 « 27r/n, the solid angle that 
the tetrahedron subtends, can be found by making use 
of L'Huilier's theorem: 



An 
7s' 



The area of each facet is: 
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where R is the radius of the sphere that covers the as- 
teroid. The height of the tetrahedron is: 



h^R(l--) 



(6) 



The volume that is enclosed between the facet and the 
sphere is: 

R 3 n - s ■ h _ V3i? 3 4 _ R 3 n 2 

3 48 _ 3V3' ( } 

The total difference in the volume between the asteroid 
and the sphere is roughly 2n ■ V: 

4tt , 4tt , 8?r 2 i? 3 

The ratio between the radii is: 

T 1 ^ 7r (9) 
R 3^/3n 

By equating the volume of the asteroid with the volume 
of a sphere with radius r we have: 



/i ■A _ 47rr 3 



and eq.© and eq.((6|) yield: 



3=1 



Therefore, we obtain: 



d s cx — . 

n 



(10) 



(11) 



(12) 



The angle, 7, between the normal of a non-equilateral 
facet and the vector joining the sphere's center to the 
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facet's centroid, is of order 9, so: 

7 cx-^. (13) 

The torque that is produced by a single facet is: 

Tj oc Sj'Yj oc n~i. (14) 

Naively, we might expect that the total torque is a ran- 
dom walk summation over the individual facets, and 
therefore 

T navve = y/nTj (X -. (15) 

n 

However, the torques from neighboring facets are not in- 
dependent. This dependence leads to a partial cancella- 
tion of the torque produced by neighboring facets that 
share a point. The torque is given by 

Toc(a + b + c)x (bxa + cxb + axc)- cos A (16) 

where A is the angle between the normal to the facet and 
the sun's flux, a, b and c are vectors from the origin to 
the vertexes located on the surface of the sphere. Con- 
sider a displacement of a single vertex, either radially or 
tangentially on the surface of the sphere, by an amount 
that would create a change of order unity in the torque, 
then the difference between the initial torque and the 
torque after displacing a by an amount da is: 

dax (bxa+cxb + axc) - cos A 

+ (a + b + c) x ((b — c) x da) • cos A 

+ (a + b + c) x (bxa + cxb + axc) - d(cos A). 

(17) 

If we displace a by an angle |da| ~ n~ 5 , then the first two 
terms in eq. (fT?| scale like 0(n~^) while the third term 
scales like 0(n~ 2 ). Likewise, if we displace a radially by 
|da| ~ n^ 1 the second term will scale like 0{n~^) while 
the first and third term will scale like 0(n~ 2 ). 

Both the radial and tangential displacement of the 
vertex a caused a change in the torque of order unity 
(C(n~2 )). However, the same vertex a is shared by sev- 
eral neighboring tetrahedrons. Summing the changes in 
the torque over all of the facets which share the vertex 
a, we find that all contributions up to 0(n~^) cancel 
out and we are left with a contribution of (D(n~ 2 ). The 
effective contribution of each such group will therefore 
scale like 0(n~ 2 ). A random walk process will now give: 

rocdf. (18) 

Following (|Goldreich fc Sarill2009l GS hereafter), we de- 
fine fy to be the ratio between the actual torque and the 
torque that would be exerted if all of the received radia- 
tion were emitted tangentially from the body's equator: 

T = "^KRHfy (19) 

where the 2/3 arises from assuming isotropic emission. 
We therefore find: 

fy~df. (20) 

We have shown that the total YORP effect is com- 
parable to the effect of a single facet. Therefore, the 



uncertainties in the asteroid's shape will induce errors in 
the torque estimate by an amount of: 

(21) 

t (R/N) d s R { ' 

where Ai? is the physical length scale of the modeling 
error and R is the length scale of the asteroid. 

4. DETAILED CALCULATION & THERMAL LAG 

In order to compute the total average torque, we first 
compute the average torque arising from each facet's 
emission independently. Let n be the unit vector point- 
ing from the facet to the sun, ds the vector perpendicular 
to the facet with a magnitude equal to its area and r the 
vector pointing from the asteroid's CM to the centroid 
of the facet, then the torque of each facet due to Lam- 
bertian reflection while neglecting specular reflection is: 

2 A<Z>d 2 , 

T re /; ee t=-- D2 (n-ds)(rxds) (22) 

where D is the distance between the sun and the aster- 
oid, A is the asteroid's albedo that is assumed to be con- 
stant. Thermal lag due to a finite thermal conductivity, k, 
might influence the obliquity change rate (it can also re- 
verse sign in extreme cases) but it hardly effects the spin 
change rate (jCapek fc Vokrouhlickvll2004h . The effect of 
thermal lag can be estimated by having some constant 
temperature per facet, T eq , and a time varying compo- 
nent, AT, which lags in time relative to the solar inso- 
lation. A facet with a constant temperature produces a 
torque which is a constant in the asteroid frame. Aver- 
aging of this torque over the asteroid's revolution around 
itself leaves us only with the projection of the torque on 
the spin axis. Therefore, the contribution of the con- 
stant temperature to the obliquity term vanishes. The 
spin component has no preference to the phase at which 
the emission takes place, hence the spin component is 
unaffected by thermal conductivity. 

As the thermal conductivity i ncreases, the ratio 
T eq /AT increases (Rubincaml I1995D . For a typical as- 
teroid, if k ■ s > 6 • 10 -5 Wj (m s K) then the equi- 
librium temperature dominates over the time varying 
component. Recent studies have shown that some as- 
teroids could hav e high enough thermal conductivity 
(jOpeil et alJl2010t l. For these asteroids, one can neglect 
the obliquity term due to thermal re-emission. In this 
paper we explore two extremes which bound the possi- 
ble behaviors, the first is the Rubincam approximation 
(k = 0) and the second is the high thermal conductivity 
regime. 

The torque due to thermal re-emission is: 

2(1- A)$d? , 

remission = — =^ ii • ds r x ds . (23) 

3 

For the high thermal conductivity regime we evaluate 
only the spin component of this torque. 

The general case of non-zero obliquity requires averag- 
ing the torque over both: 

(i) The orbit of the asteroid around the sun. 

(ii) A revolution of the asteroid around itself. 
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For a zero obliquity orbit, it is sufficient to average over 
either one of the above. Since we assumed that the spin 
period and the orbital period are non-commensurate, we 
can calculate these averages in arbitrary order. We find 
that it is more convenient to first average over the orbit 
around the sun, while holding each facet pointing in a 
fixed direction. 

4.1. Heliocentric Orbit Average 

The insolation, averaged o ver the orbi t of the asteroid 
around the sun, is given by (|Wardlll974D 



and, 



1 

< I >= - 



L(7) ~ $sin(0) . , 
n-dsdt= — (24) 



T J AttcD 2 

where T is the heliocentric orbital period and is the 
angle between ds and the normal to the orbital plane. 

4.2. Asteroid Spin Averaging 

Since we are interested in calculating the torque in 
the inertial system (system 1), we need to transform the 
torque from the asteroid system (system 3) to the iner- 
tial system. In the asteroid system the torque per facet 
is: 



/ / 

T = T 

3 



reflection^ 



emission, j 



V x < J > ds' 

3' 



2S, 



y'j cos 9 j — z'j sin tTj sin ^ 3 
< I > | z'j sin 9, cos - x[. cos 9'. 

, sin 0' (x\j sin (pj - Vj cos (p, ) 



(25) 



where Sj is the facet's area, and 9'j, <p'j are the polar and 
azimuthal angles accordingly. The torque is then trans- 
formed into system 1 with the standard Euler angle ro- 
tation matrix: 

C— cos e sin -0 — cose cos ip sine\ 
cos^f -sinip (26) 

sin e sin tp sin e cos tp cos e / 

where tp is the rotation angle of the asteroid around itself 
and e is the obliquity. For tp = the asteroid system's 
x' axis coincides with the y axis of the inertial system. 
The insolation in system 1 is: 



$sin(0j) $sin(cos _1 ([R • ds] • z)) 



7T 
' \2 



_ $[(cos(V> + (/>' j ) sin O'j) 
+(cos0' sine — cosesin^ + (pj) sin0') 2 ]5 



(27) 



7T 

The averaged torque in the inertial system can now be 
rewritten as: 

OTT Jo 

In order to calculate the torque we define the following 
numerical functions: 



B(9',e)= — J [(cos V sin 9'f 



(29) 



+ (cos 0' sin e — cos e sin tp sin 0') 2 ] 2 dtp 



1 C 2lx 

G(6 , ,e) = —s sin^[(cosV'sin0') 2 
K ' 2ir 2 Jo (30) 

-I- (cos 0' sine — cos e sin i/; sin 0') 2 ] ?dip. 

Notice that the physical interpretation of the function B 
is simply the average insolation a facet receives at a given 
angle 9' . T his result is similar to the one deriv ed in lWardl 
(| 19 741) and iNesvornv k. Vokrouhlickvl (|2007l ). Figured 
shows B and G as function of 0' for various values of e. 

Rewriting eg. (|2"51) with eq. (|2"9"|) . ea. (j3U]) and using 
eq. (f26j) and eq. (|25|) . we obtain by summing over all of 
the facets ds,, (assuming n facets): 



3 

3=1 



{B(9'j , e) sin e sin 9'j (x'j sin <p'j — y'j cos <pL ) 



+ G(9'j,e) cos e cos 0^- (:Ej sin <pj — y'j cos ^ ) } 
n 2 9 <i> 

Ty = -J2 -^-G(d' j ,e){co S d , j (x' j cos $ + y'j sin ^. ) 



3=1 

z'sin0^} 



™ 2 S 1 

7i = ~^-{B(0^e) cosesin0;( a; ; sin^ - y'j cos$) 



3=1 

+ G(0j , e) sin e cos 9'j {—x'j sin + y'j cos <^-)} 

(31) 

We have so far treated the asteroid as having homo- 
geneous density, so that the center of mass is at the ge- 
ometric center. Inhomogeneous density would result in 
a displacement of the center of mass from the geometric 
center and can easily be treated by adding the displace- 
ment term, a vector pointing from the geometric CM 
toward the true CM, to r' in eq.t|25|) . 



4.3. Spin Evolution 
In the inertial system the unit spin vector is: 




(32) 



The rate of change of the spin vector's magnitude is: 



C 



(33) 



where C is the asteroid's moment of inertia. 
Substituting eq. (|3"Tj) into eq. (|3"3")) we obtain: 

" 2S' J $£(0',e)sin0'. , , 
< s >= - 2^ ^ -Oj sin ^ " Vj cos (pj ) . 



3=1 



(34) 

This result depends only on B since the torque along 
the z' axis is the only vector that is fixed in all of the 
coordinate systems. It has the simple physical interpre- 
tation as the average insolation times the torque that 
each facet contributes. Thermal conductivity has no 
effect on the spin evolution. An interesting result is 
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Fig. 1. — The functions B and G for different obliquities as a 
function of the facet's polar angle relative to the spin axis. The 
dotted red line represents e ?s 55° . 



that when the obliquity is approximately e ~ 55°, B 
is nearly a constant with respect to 9', and as a result, 
s tends to vanish due to the same argument that was 
used to show that s vanishes due to absorption (refer to 
[Rub incam fc PaddacH (]2010f )). A clue as to why B is 
roughly 9' independent for e 55° might be found in 
the fact that a related expression: 

(cos ip sin 9') 2 + (cos 9' sine — cosesin^sin^') 2 ^ 

(35) 

is independent of 9' for cos2e = —1/3. This 
result coincides with t he v alue calculated by 
iNesvornv fc Vokrouhlickyl (|2007l ). Note however, 
the critical obliquity, e 55°, is by no means a stable 
point. Even though s = 0, the obliquity will evolve away 
from the critical obliquity (see £|4.4[) . Our calculations 
were tested by computing the YORP effect on Kleopatra 
and 1998KY26, presented in FigJH and comparing it 
with the results of ICapek fc Vokrouhlickyl (|2004[ ) in 
Fig. 3 of their paper for 1998KY26. Our results are 
about 5% smaller than theirs. 

4.4. Obliquity Evolution 
The obliquity change rate can be derived as follows. 



From eg . ([32]) we have: 




1 1.5 2 

Obliquity(rad) 

Fig. 2.— The YORP effect calculated for 1998KY26 and Kleopa- 
tra. The solid line represents the spin evolution coefficient fy and 
the dashed line represents the obliquity evolution coefficient nor- 
malized in the same fashion as fy. 1998KY26 is a typical Type I 
asteroid while Kleopatra is a typical Type II asteroid. Zero thermal 
lag is assumed. 

Differentiating with respect to time we obtain: 

s • z (s • z)s 



and rearranging yields: 

t\s cos e — z] 



t x cos e — t z sin e 



Cs sin e 

Substituting in eq.(|23]l we get: 

" 25 J $G(6»',e) cos &, 



Cs 



(37) 



(38) 



3Cs 



cos e 



s • z 

s 



(36) 



(x'j sin cji'j — y'j cos 4>'j)Q 
(39) 

where = 1 for the Rubincam approximation and O = 
A for the high thermal conductivity regime. Just like the 
function B governs the spin's evolution, G governs the 
obliquity's evolution. 

4.5. Comparison of Qualitative with Precise Results 

In order to check our scaling results of section 2500 
random bodies were constructed in the manner described 
in <|3] For each body we computed fy and d s . Figure [3] 
shows fy for the randomly constructed bodies as a func- 
tion of spherical deviation, d s , along with the best fit for 
cq. (|20|) . In addition, fy for 18 real asteroids (their shape 
was taken from ECHO JPL) was computed using our al- 
gorithm and plotted together with the 2 measured YORP 
for comparison. For most of the asteroids we can see a 
good correlation between our qualitative results and our 
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Fig. 3. — Maximal /y- vs. deviation from sphere, the red line 

3 

denotes the best fit /y = 0.14 • d£ for the simulated asteroids 
and the dashed red line denotes the best fit for the real asteroids 

3 

fy = 0.04dg ; zero eccentricity is assumed. The filled red stars 
represent real asteroids calculated with eg. (I34t . The filled black 
stars are the computed /y for measured asteroids (normalized to 
take into account their current obliquity and eccentricity) and the 
hollow circles are their measured YORP respectively. The asteroids 
that were calculated are: 216 Kleopatra. 1620 Geographos, 1992 
SK, 1998 KY26, 52670 1998 ML14, 2002 CE26, 2053 Bacchus, 
4179 Toutatis, 4486 Mithra, 4769 Castalia, 6489 Golveka, 25143 
Itokawa, 433 Eros, 1580 Betulia, 54509 YORP (2000 PH5), 1999 
KW4 (alpha+beta) and 2100 Rashalom. Zero thermal lag was 
assumed. 

precise calculations, although our random asteroids tend 
to have a higher fy. The ratio between our calculated 
results and the measured YORP for 1620 Geographos 
and 2000 PH5 are 2.18 and 5.22 respectively. Taking the 
mean of the real asteroids, we find that fy ~ 2.8 ■ 10~ 3 . 
Reexamining the order of magnitude approach that was 
taken by GS, we can see that their choice of fy 6 x 10~ 4 
is a bit low but is still within acceptable range. 

5. EVOLUTION 

Since under the mapping z to —z, the obliquity chang- 
ing component of the torque does not change sign but 
the z component of the spin vector does, we have the 
antisymmetry: G(8,e) = — G(8, ir — e). Similarly, since 
the spin changing component of the torque does change 
sign: B(8,e) = B(6,n — e). It can be seen from eq. (|3l?f 
that every asteroid, regardless of its shape, has at least 
two equilibrium points with regard to the change in the 
obliquity since for e = or tt/2, G is zero. Whether they 
are stable or not depends on the derivative of ea. (l3T)l) 
with respect to e, which depends on the asteroid's ge- 
ometry. An asteroid will have an identical number of 
equilibrium points in the range [0,7r/2) and the range 
(7r/2,7r] due to the symmetry of G. The stability of e at 
e = determines the stability of e at all the other equilib- 
rium points since the equilibrium points change stability 
alternatingly. iVokrouhlickv fc CapekJ (|2002f) have classi- 
fied the two most frequent spin and obliquity behaviors 
of asteroids into Type I and Type II. Type I is catego- 
rized by having e > for < e < tt/2 and a spin change 
rate that is positive for e = and negative for e = 7r/2. 
Type II is the opposite of Type I, categorized by hav- 
ing e<0for0<e<7r/2 and a spin change rate that 
is negative for e = and positive for e = tt/2. Both 
Type I and II are defined to have no nodes for e in the 
range (0, 7r/2). Our calculations show that 78% of our 



randomly drawn asteroids are divided into Type I and II 
with equal likelihood. Our calculations also show that 5 
out of the 18 real asteroids are not Type I or II. Of the 
remaining 13 asteroids, 9 were Type I and 4 were Type 
II, which is consistent with our random asteroids. One 
of the curious aspects of both types, is that their spin 
tends to slow down at their stable point with respect to 
obliquity. An asteroid which starts with an obliquity in 
the range of (0, tt/2), will tend to evolve toward e = 
if it is a Type II asteroid where its spin change rate is 
negative. A Type I asteroid will evolve toward e = it/ 2 
where its spin change rate is also negative. Notice that 
once the spin rate reaches zero, the asteroid does not sim- 
ply change its obliquity from e to 7r — e since that would 
also involve inverting the z' axis. Rather the obliquity 
remains fixed while G changes its sign. This causes the 
equilibrium point to lose its stability and results in the 
asteroid evolving to its other equilibrium point where the 
spin will once again slow down. 

A viable explanation for this phenomena can be found 
by inspecting the stability at e = 0: 



de 
de 



e=0 



E 



3irCs 



cos 8j ( Sj x'j sin 0' — Sj y' cos 0' ) . 



The spin change rate with zero obliquity is: 
2$ 



(40) 



s(e = 0) = - 



^ 3irC 



sin' 



6j ( Sj x'j sin (j>'j — Sj y '• cos <j>'j ) . 

(41) 

The correlation can be understood as follows. In $3] 
we saw that there is a strong anti-correlation between 
neighboring facets that causes a reduction in the YORP 
magnitude. This anti-correlation is represented by the 
expression in the parentheses in eq.(j40|) and eq. (j4Tj) . To 
mimic this anti-correlation we take a set, {i n }, of N num- 
bers randomly drawn uniformly from -1 to 1 and repre- 
sent the term in the parentheses as the difference between 
any two consecutive numbers. For example for the spin 
change rate we can write: 



JV-l 

s(e = Q) = - (x n+1 - x n ) sin 2 8' n , 

n=l 

and for the obliquity stability: 

\ (.t„ + i - x n ) cos 2 e' n 



de 
Je 



N-l 

E 1 

n=l 



(42) 



(43) 



where 8' is a randomly drawn vector of length N — 1 
ranging from to n and sorted in an increasing manner. 
Since in asteroids we have 



Sj sin 0' (x' sin $ - y\ cos <f{. ) = 0, (44) 



we set x N - x N -i = Y, n =i ( x n+i - x n )sin8' n /sind N ^i. 

The sign of <ie/rfe| e= o ■ s(e = 0) for our simple model 
has a likelihood of 85% being positive, while the likeli- 
hood of de/de\ e =o ■ s(e = 0) for our randomly constructed 
asteroids being positive is 88% (it is more than 78% be- 
cause it accounts for cases with more than only 2 stable 
obliquity points). 
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The strong tendency of asteroids to spin-down is 
merely a result of the correlations between neighboring 
facets. 

Asteroids with moderate thermal conductivity might 
exhibit different obliquity behavior which is beyond the 
scope of this work (jCapek fc Vokrouhlickvl[200l . 

When the spin of the asteroid is sufficiently low, the 
asteroid might fall into a chaotic tumbling rotation state 
and it is not clear if the asteroid can recover principal 
axes rotation (|Vokrouhlickv et al]l2007h . 

If the asteroid is a rubble pile, the spin-up might even- 
tually break up the asteroid due to centrifugal force. 
However, ru bble-pile asteroid deformations are not well 
unde rstood (|Holsappldl20lol IScheeresll200l iWalsh et all 
120081). Inspectio n of the shape of binary 1999KW4 
(jOstro et al.ll2006T ) suggests that break up is possible. 

6. BYORP 

The BYORP effect is very similar to the YORP effect 
so the understanding of YORP continues into BYORP. 
We will assume that the secondary is synchronously 
locked with the primary, and that the secondary's spin 
vector is aligned with the orbital angular momentum. 
This assumption is reasonable due to the strong tidal in- 
teractions between the binary components. The primary 
is assumed not to be locked with the secondary (i.e., not 
a double synchronous state), and thus the evolution of 
the binary orbit does not depend on the orientation of 
the primary. In the following discussion, the origin of the 
asteroid system (system 3) will be at the CM of the sec- 
ondary asteroid with its smallest principle axis, labeled 
the x' axis, pointing away from the primary. The iner- 
tial system's (system 1) origin will now be located at the 
CM of the primary, the z axis aligned with the orbital 
angular momentum and the x axis pointing toward the 
periapsis of the binary's orbit. For the special case where 
the orbit is circular, we define the x axis to point along 
the line of the ascending node. The angle ip is defined 
relative to the line of nodes where our plane of reference 
is the orbit of the primary around the sun. In order to 
be consistent with our previous definitions, the rotation 
angle ip is now defined as: ip = A + w, where A is the 
mean anomaly and w is the argument of periapsis. In 
addition, the role of the inclination of the binary orbit 
relative to the primary's orbit around the sun in BYORP 
takes the role of the obliquity in YORP, e f-> i. 

In YORP the main interest is to find the evolution of 
the obliquity and spin. In BYORP we are interested in 
the evolution of the binary orbit's inclination relative to 
the orbit of the primary around the sun, the evolution 
of the binary semi-major axis and the evolution of the 
binary's eccentricity. Since the secondary is assumed to 
be locked, changes in the orbital spin vector correspond 
to change in the semi-major axis, while the change in 
the inclination corresponds to the change in obliquity in 
YORP. There is no analog for the eccentricity change in 
YORP. 

Our derivation of BYORP neglects the distance of the 
secondary's facets from its CM and we assume that all of 
the facets are located at the CM of the secondary. This 
is justified since in typical cases the error in the torque 
due to errors in the facets' locations will be at least an 
order of magnitude less than the actual torque. 

Just like in YORP, we will inspect the same two ex- 



treme of very high and very low thermal conductivity. In 
the following section the eccentricity of the orbit of the 
secondary around the primary will be denoted e, to dis- 
tinguish from the eccentricity of the binary orbit around 
the sun e. 

Absorbtion of radiation leads to no secular changes in 
the semi-major axis, eccentricity and inclination of the 
binary. This can be shown us ing a similar argument as in 
iRubincam k, Paddackl(|2010D . Fixing the binary's orbital 
elements, the torque at a given point along the heliocen- 
tric orbit cancels out with the torque arising at a point 
with the true anomaly of the heliocentric orbit changed 
by 180°, since the radiation force changes sign. This is 
true for any point along the orbit, so the net effect of 
absorption vanishes. 

6.1. Semi-Major Axis Evolution 

The change in the semi-major axis can easily be de- 
duced from the change in the energy of a binary orbit. 
To zero order in the eccentricity the change is: 



where a is the binary's semi-major axis , n is the mean 
orbital angular velocity and Ft is the transverse orbital 
force per unit mass arising due to BYORP, positive in 
the direction of movement. In order to compute < a > 
we define two rotation matrices. The first, Rf , describes 
a rotation around the z axis in the inertial frame by an 
angle / + w, where / is the true anomaly. The second 
rotation matrix, Ra, describes a rotation around the z 
axis by an angle ip. With our choice of coordinates: 

n 9 q . 

F > = - E ut- < 1 > ^ ■ ds V ■ ( R f • y') ( 46 ) 

here M s is the secondary's mass. By first averaging over 
the heliocentric orbit and then averaging over the binary 
orbit, we obtain: 

< F * >= ~ E 3 a/ sin 6 3 sin #i ■ ( 47 ) 

i=i s 

This result is intuitive since the expression merely sums 
the transverse force from each facet. Eq. (|47|) resembles 
the YORP spin change rate in eq. ([34]) . Just like the spin 
change rate for YORP was unaffected by thermal lag, 
the semi-major axis change rate is unaffected by thermal 
lag as well. When i rj 55°, < d > vanishes for the same 
reasons that < s > vanishes. We define /by, similarly 
to fy, to be the ratio between the BYORP actual effect 
and the effect it would have on the semi-major change 
rate if all of the received radiation had been emitted tan- 
gentially from the secondary. Therefore, the torque can 
be rewritten as: 

t = \aR 2 f BY <S>. (48) 

In the derivation in section fJ3]the angle between the force 

i_ 

and the lever that produced the torque scaled like 0(di). 
However, in BYORP there is no correlation between the 
force and the lever so we would expect: 

fBY-d s ^f Y . (49) 




Fig. 4. — The ratio of fgy to fy- The central line corresponds 

2 

to the best fit of fsY = 0.43/y , while the long dot dashed line 

2 

represents the upper bound of Jby = 1-85/y and the short dashed 

2 

line represents the lower bound of /by = 0.1 f Y of the 2a disper- 
sion. Also plotted is the ratio for our real asteroid models, using 
the same asteroids listed in Fig|3] Zero thermal lag is assumed 

Fitting the ratio between the coefficients with /by = 

2 

c- fy we find that the ratio for our randomly constructed 
bodies satisfies: 



f BY = 0.43/|. 



(50) 

The 2(7 dispersion around the best fit number has an 

2 

upper limit of /by = 1-85/y- and a lower limit of Jby 

■ 

0.1/y. Figure S] shows the ratio for our random asteroids 
as well as for our real asteroid models. Taking the mean 
of the real asteroid's Jby yields: 

f BY « 0.01. (51) 

The uncertainties in the asteroid's shape will induce 
errors in the torque estimate by an amount of: 



Ax 

T 



2AR 



2AR 



3(R S /N) 3d s R s 



(52) 



where AR is the physical length scale of the modeling 
error, R s is the length scale of the secondary and the 
factor 2/3 arises from the relation of YORP to BYORP 
strength. 

6.2. Inclination Evolution 

6.2.1. Zero Order 

The e volution of th e inclination can be calculated from 
(see e.g. IBurnslH976f) . 

4- = ku _1 (l -e 2 )]^F N cos{f + w)/{\ + ecos/) (53) 
at 

where Fn is the normal component of the BYORP force 
per unit mass, relative to the binary orbit. 



Fn — — 



v -v 2Sj 

3 = 1 



■ < I > (Ra • ds',) ■ (Rf • z 1 ) (54) 



Due to the quadrapole moment of the primary, char- 
acterized by its J 2 , the angular momentum of the sec- 
ondary will be coupled to the angular momentum of the 
primary ()Cuk fc Nesvorny|2010D . This interaction causes 



the orbit of the secondary to rapidly precess, typically at 
timescales of a few months. This coupling requires us 
to modify the inclination's evolution rate by a factor j3, 
the ratio of the secondary's orbital angular momentum 
relative to the total angular momentum (which may be 
dominated by the primary's spin). In order to take into 
account the rapid precession, an averaging of the argu- 
ment of periapsis is required. No averaging of the line of 
ascending node is necessary since it does not enter our 
equations. 

The sun will also cause the binary orbit to precess. 
However, since the timescale for this interaction is of the 
order of hundreds of years, the angular momentum of the 
secondary's orbit will stay aligned with the primary's 
angular momentum. This precession will introduce no 
additional change in our equations since heliocentric line 
of nodes does not enter our equations. 

The average inclination change rate, up to zero order 
with respect to the eccentricity is: 



ali 

< dt > " 



^2S J G{9' J ,i)<S> 



3=1 



3Af,na 



9/3 cos 9'j sin fy 



(55) 



where we added the term just like in the case for the 
obliquity evolution. This result resembles eq. ([39l) . 

6.2.2. First Order 

In order to simplify our calculations we expand / to 
first order in the eccentricity, / A + 2e'sinA. To 0(e) 
order in eccentricity the contribution is: 



ali Si e$ 
< — > = > — 

dt 3M s na 

3=1 



cos 0' {35(0', i) cos iu- QHieLi) 



x cos(2</(/. + w)} 



(56) 



(57) 



where the function H is defined to be: 

i r 27T 

H(9,i) = — - / [cos 2V»(cos?/>sin0) 2 
^ Jo 

+ (cos0sini — cosisin^sh^) 2 ]^?/; 

Averaging over the precession of the argument of peri- 
apsis cancels the first order correction due to eccentricity. 

The reasoning presented in pleads us to expect a cor- 
relation between the change in the inclination and the 
change in the semi-major axis. 

6.3. Eccentricity Evolution 

Unlike the change in the semi-major axis and inclina- 
tion which have their analogies in YORP, there is no 
analogy for the eccentricity evolution. The eccentricity 
vector is defined as: 



v x h 



(58) 



where h is the binary's orbital angular momentum per 
unit mass, /x = GM p and v is the orbital velocity of the 
secondary. Since in an unperturbed orbit this vector is 
conserved, its time derivative is: 



^ BYORP X h V X T BYORP 



/' 



(59) 
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where Fbyorp is the force that arises from BYORP and 
tbyorp is the torque that is produced by BYORP. The 
radius vector r' in the asteroid system is: 



(60) 



The time derivative of the eccentricity vector can be writ- 
ten as: 




_ n o<f) q 

J2=i^<I>{(K x ds)xh 



3=1 

d(R f r') 

dt 



3// 



(61) 



x ((R f r') x (R x ds ))}. 



Averaging this result over the binary's period, expanding 
up to O(e') and taking the projection on the secondary's 
orbital plane yields: 



< e T , > = 



2, S^sinfl' f 

3 \ 2QG(6' j , i) [cos(2^- + w) - 3 cos w) 



E 



6M s na 

3 = 1 

+ e[2B(e' j ,i) smc^'j - eH(9' j , i)(sin(<^ + 2w) 



+ 3sin(3^. +2w))] 

2eG(^,z)[sin(20;.+ W ) 
bM s na J J 



" sin 0'. 
E 

3=1 



- 3 sin w)} + e[-AB(e'j,i) cos (p'j 

+ QH(9' p i)(3 008(3$ + 2w) + cos($ + 2w))} 

(62) 



6.3.1. Eccentricity's Magnitude Evolution 

The change in the magnitude of the eccentricity is the 
projection of e on the x axis. To zero order in eccentricity 
this change is: 

" SiG(6Lm sin 6>' r , , 

< e >= V 3 3 -9 cos(2<£' + w) - 3 cosw] 

^— ' 6M s na J 



3=1 



This contribution is zero for i = 0,7r/2. 
Up to 0(e) the change is: 



(63) 



< e > 



3 -e [2B (d'^i) sin0 - 0^,* 

bM s na j j j 



E 

3=1 

x (sin(0' + 2w) + 3 sin(30' + 2w))]. 



(64) 

Averaging over the precession of the argument of periap- 
sis retains only the term: 



< e >= E aS^ 58 ^' sin ^ sin ^- (65) 

3 = 1 



Compaxing eq. (|45|) and eq. (|47)) to eq.([65|) we find, simi- 
larly to GS, that the eccentricity evolves like: 



o.oi 



0.005 



J o 



-0.005 



-0.01 



* \ 

* \ 
t \ 
t \ 

t 

1 


\ / 
\ / 
\ / 
\ / 
% / 

\ / 

\ \ / * 
\ \ / / 
\ x / * 
\ / * 

\ / * * 
\. S s * 


di/dt 




da /dt 





0.5 



1 1.5 2 
Inclination(rad) 



2.5 



Fig. 5.— The BYORP effect calculated for 1999KW4. The solid 
line is the semi-major axis changing coefficient Jby and the dashed 
line represents the inclination changing torque coefficient normal- 
ized in the same fashion as Jby- 1999KW4 is a typical Type II 
asteroid, zero thermal lag is assumed. 



e (X a 4 . 



(66) 



We recover the result of iMcMahon fc Scheerei (|2010f ) 
that circular orbits remain circular. An interesting result 
is that both the semi-major axis evolution and eccentric- 
ity evolution are independent of the asteroid's albedo and 
thermal conductivity, in the two extreme regimes of high 
and low thermal conductivities. 

6.4. Results 

Our calculations were tested on the binary asteroid 
1999KW4 which is the best modeled binary. Figure 
[5] shows our calculated change in the semi-major axis 
and in the inclination of the binary. Our calculations 
show that the secondary is currently drifting away from 
the primary at a rate of about a = 7 cm/year, which 
is in an agreement to wit hin 5% with the result of 
IMcMahon & Scheeresl pOloh . 

Due to the similarity between YORP and BYORP, the 
statistics found in Sj5] are similar for BYORP. We pre- 
formed calculations on the 2500 randomly constructed 
asteroids. The stability of di(i = 0)/dt with respect to a 
change in the inclination determines the stability of the 
rest of the equilibrium points. The correlation that was 
found between the sign of s(e — 0) and the stability at 
e = also holds for the sign of a(i = 0) and the sta- 
bility at i = 0. For randomly shaped asteroids, there is 
an equal likelihood for i — to be a stable point or an 
unstable point and 73% were either Type I or Type II 
(with equal likelihood). Our calculations also show that 
2 out of the 18 real asteroids are not Type I or II. Of the 
remaining 16 asteroids, 6 were Type I and 10 were Type 
II, which is consistent with our random asteroids. 

Unlike YORP, the timescale for reaching the stable 
point can be longer than the lifetime of the system, so 
not every binary will reach it. In order for the system to 
reach its stable point, the difference between the starting 
inclination and the inclination at the stable point needs 
to be Ai « 

7. CONCLUSIONS 

In this paper we have developed an intuitive and sim- 
ple analytical model that predicts the behavior of the 
YORP and BYORP effects. Our model was computed 
on randomly constructed asteroids and on shape models 
of real asteroids. The calculations of our simple model 
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show that the YORP effect is dominated by the inner cor- 
relations between facets. These correlations explain the 
magnitude of the effect and its preference to spin-down 
the asteroid. We have also shown that the dimensionless 
parameter of d s governs the strength of the effect (but 
not of the sign). 

Kilometer sized or smaller near earth asteroids undergo 
a relative fast spin-up or spin-down due to YORP. The 
YORP effect tends, in most cases, to change the initial 
obliquity of the asteroid to cither e = 0ore = 7r/2. For 
both extreme regimes of thermal conductivity asteroids 
will most likely spin-down. The rapid spin-up of aster- 
oids, will eventually cause the asteroid to break up and 
possibly form a binary. The binary will evolve under 
BYORP and tidal forces; the former will tend to orient 



the binary into a i — 0ori~7r/2 orbit, if the system is 
near enough relative to the stable point. After the binary 
orbit reaches its stable inclination, the binary will most 
likely migrate inwards until it is stopped by tidal forces 
(if it survives the migration). A simple relation between 
the magnitude of the YORP effect and BYORP effect 
is calculated, thus knowledge of the former can help us 
estimate the strength of the latter. 

We would like to thank Lance Benner for providing us 
the shape models for the real asteroids. We would also 
like to thank the referee Matija Cuk for alerting our 
attention about the importance of the J2 interactions. 
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